%load_ext autoreload
#load in modules
%autoreload 2
from wildfireassessment.ops import * #my package
import rasterio
from rasterio.transform import xy
from rasterio.plot import show, show_hist
from rasterio.mask import mask
import matplotlib.pyplot as plt
from pathlib import Path
import os
import earthpy as et
import earthpy.plot as ep
import numpy as np
import geopandas as gpd
from fiona.crs import from_epsg
from pyproj import Proj, transform
from functools import partial
import pyproj
from shapely.ops import transform
from rasterstats import zonal_stats
import itertools
%matplotlib inline
# from osgeo import gdal, osr, ogr, gdal_array
# from osgeo.gdalconst import *
#read in filepaths for data
filepath_post = Path("./data/Paradise/post")
filepath_pre = Path("./data/Paradise/pre")
#WorldView Post/Pre
fps_wv_post = list(filepath_post.glob("2*_clip.tif"))
fps_wv_pre = list(filepath_pre.glob("2*_clip.tif"))
#Sent2 Post/Pre
fp_sent2_post = filepath_post / "B08_post_clipped_byte_scaled.tif"
fp_sent2_pre = filepath_pre / "B08_pre_clipped_byte_scaled.tif"
raster_src_post, rgb_post = readRGBImg(fps_wv_post[5])
raster_src_pre, rgb_pre = readRGBImg(fps_wv_pre[5])
raster_src_post_b08, b08_post = readOneImg(fp_sent2_post)
raster_src_pre_b08, b08_pre = readOneImg(fp_sent2_pre)
!gdalinfo data/Paradise/post/B08_post_clipped_byte_scaled.tif
fig, (ax1, ax2) = plt.subplots(ncols=2, figsize=(15, 10))
im = ax1.imshow(rgb_pre)
ax1.set_title("Pre-fire")
im = ax2.imshow(rgb_post)
ax2.set_title("Post-fire")
plt.show()
x_r_pre = rgb_pre[:,:,0].ravel().reshape(-1, 1)
x_g_pre = rgb_pre[:,:,1].ravel().reshape(-1, 1)
x_b_pre = rgb_pre[:,:,2].ravel().reshape(-1, 1)
x_pre = np.hstack((x_r_pre, x_g_pre, x_b_pre))
x_r_pre, x_g_pre, x_b_pre = None, None, None
y_r_post = rgb_post[:,:,0].ravel().reshape(-1, 1)
y_g_post = rgb_post[:,:,1].ravel().reshape(-1, 1)
y_b_post = rgb_post[:,:,2].ravel().reshape(-1, 1)
y_post = np.hstack((y_r_post, y_g_post, y_b_post))
y_r_post, y_g_post, y_b_post = None, None, None
n_bins = 20
colors = ['red', 'green', 'blue']
fig, axs = plt.subplots(1, 2, sharey=True, tight_layout=True, figsize=(15, 10))
axs[0].hist(x_pre, n_bins, alpha=0.5, color=colors, label=colors)
axs[1].hist(y_post, n_bins, alpha=0.5, color=colors, label=colors)
axs[0].legend(prop={'size': 10})
axs[1].legend(prop={'size': 10})
axs[0].set_title("Pre-fire")
axs[1].set_title("Post-fire")
plt.show()
MinMax histogram stretch on values on Sent2 images. The maximum was determined in properties of pre-fire image. This is done since DigitalGlobe images were adjusted according to their Dynamic Range Adjustment (no specifics on how this histogram stretch is done).
#command in gdal translate
#!gdal_translate -scale 0 7210 0 255 -ot Byte data/Paradise/post/B08_post_clipped.tif data/Paradise/post/B08_post_clipped_byte_scaled.tif
#command in gdal translate
#!gdal_translate -scale 0 7210 0 255 -ot Byte data/Paradise/pre/B08_pre_clipped.tif data/Paradise/pre/B08_pre_clipped_byte_scaled.tif
bbox_post = box(raster_src_post.bounds.left, raster_src_post.bounds.bottom, raster_src_post.bounds.right, raster_src_post.bounds.top)
bbox_pre = box(raster_src_pre.bounds.left, raster_src_pre.bounds.bottom, raster_src_pre.bounds.right, raster_src_pre.bounds.top)
out_img_post_b08, out_img_transform_post_b08, out_meta_post = clipImg(bbox_post, raster_src_post_b08, proj=4326)
out_img_pre_b08, out_img_transform_pre_b08 , out_meta_pre = clipImg(bbox_pre, raster_src_pre_b08, proj=4326)
fig, (ax1, ax2) = plt.subplots(ncols=2, figsize=(15, 10))
im = ax1.imshow(out_img_pre_b08[0], cmap='gray')
ax1.set_title("Pre-fire")
im = ax2.imshow(out_img_post_b08[0], cmap='gray')
ax2.set_title("Post-fire")
plt.show()
x_nir_pre = out_img_pre_b08[0].ravel()
y_nir_post = out_img_post_b08[0].ravel()
n_bins = 20
colors = ['red']
fig, axs = plt.subplots(1, 2, sharey=True, tight_layout=True, figsize=(8, 5))
axs[0].hist(x_nir_pre, n_bins, alpha=0.5, color=colors, label='nir')
axs[1].hist(y_nir_post, n_bins, alpha=0.5, color=colors, label='nir')
axs[0].legend(prop={'size': 10})
axs[1].legend(prop={'size': 10})
axs[0].set_title("Pre-fire")
axs[1].set_title("Post-fire")
plt.show()
indices = chunkImageIndices(rgb_post)
indices
coord_ind = (9792-1024, 9792+1024, 9792-1024, 9792+1024)
from skimage.segmentation import slic, felzenszwalb
from skimage.segmentation import mark_boundaries
from skimage.measure import label, regionprops
segments_slic = slic(rgb_pre[coord_ind[0]:coord_ind[1], coord_ind[2]:coord_ind[3], :], n_segments=1000, compactness=15)
#segments_slic = slic(rgb_pre[indices[1][0]:indices[1][1], indices[1][2]:indices[1][3], :], n_segments=5000, compactness=10)
#segments_slic = np.add(np.ones(segments_slic.shape), segments_slic).astype(np.uint16)
print('SLIC number of segments: {}'.format(len(np.unique(segments_slic))))
#segments_fz = felzenszwalb(out_img_pre[0:chunksize[0], 0:chunksize[1], :], scale=10, sigma=0.5, min_size=20)
#print('Felz number of segments: {}'.format(len(np.unique(segments_fz))))
boundary_rgb_post = mark_boundaries(rgb_post[coord_ind[0]:coord_ind[1], coord_ind[2]:coord_ind[3], :], segments_slic, outline_color=(1, 1, 0), mode='thick')
bouncary_rgb_pre = mark_boundaries(rgb_pre[coord_ind[0]:coord_ind[1], coord_ind[2]:coord_ind[3], :], segments_slic, outline_color=(1, 1, 0), mode='thick')
# boundary_rgb_post = mark_boundaries(rgb_post[indices[1][0]:indices[1][1], indices[1][2]:indices[1][3], :], segments_slic, outline_color=(1, 1, 0), mode='thick')
# bouncary_rgb_pre = mark_boundaries(rgb_pre[indices[1][0]:indices[1][1], indices[1][2]:indices[1][3], :], segments_slic, outline_color=(1, 1, 0), mode='thick')
#boundary_rgb_post = mark_boundaries(out_img_pre[0:chunksize[0], 0:chunksize[1], :], segments_fz, outline_color=(1, 1, 0), mode='thick')
fig, (ax1, ax2) = plt.subplots(nrows=2, figsize=(15, 30))
im = ax1.imshow(bouncary_rgb_pre)
ax1.set_title("Pre-fire")
im = ax2.imshow(boundary_rgb_post)
ax2.set_title("Post-fire")
plt.show()
from affine import Affine
(indices[0][1]*indices[0][3]//(2048*2048))*1000
indices[1]
raster_src_pre.transform
raster_src_pre.index(raster_src_pre.transform[2], raster_src_pre.transform[5])
Affine.translation(1.0, 1.0) * raster_src_pre.transform
raster_src_pre.height
tup = raster_src_pre.xy(8168, 0)
print(tup)
tup_origin = raster_src_pre.xy(0,0)
print(tup_origin)
abs(tup[0] - tup_origin[0])
abs(tup[1] - tup_origin[1])
x_trans = abs(tup[0] - tup_origin[0])
y_trans = abs(tup[1] - tup_origin[1])
Affine.translation(x_trans, -y_trans) * raster_src_pre.transform
raster_src_post.crs
import shapely
from shapely.wkt import loads
from shapely.geometry import Polygon, mapping, shape, MultiPolygon, box
from rasterio import features
import geopandas as gpd
import pandas as pd
shapes_list = [{'seg_index': int(v), 'geometry': loads(shape(g).wkt)} for g, v in features.shapes(segments_slic.astype(np.uint16), mask=None, transform=Affine.translation(x_trans, -y_trans) * raster_src_pre.transform)]
gdf = gpd.GeoDataFrame(shapes_list)
gdf.crs = {'init': 'EPSG:4326'}
gdf.insert(2, 'crs', str(raster_src_post.crs))
gdf.head()
gdf.to_file("segments.shp")
## use regionprops onn segments to extract properties for dataframe
#post
regions_red = regionprops(segments_slic, rgb_post[0:chunksize[0], 0:chunksize[1],0])
regions_blue = regionprops(segments_slic, rgb_post[0:chunksize[0], 0:chunksize[1],1])
regions_green = regionprops(segments_slic, rgb_post[0:chunksize[0], 0:chunksize[1],2])
#pre
regions_red_pre = regionprops(segments_slic, rgb_pre[0:chunksize[0], 0:chunksize[1],0])
regions_blue_pre = regionprops(segments_slic, rgb_pre[0:chunksize[0], 0:chunksize[1],1])
regions_green_pre = regionprops(segments_slic, rgb_pre[0:chunksize[0], 0:chunksize[1],2])
gdf = gdf[gdf['seg_index'] != 0]
project = partial(
pyproj.transform,
pyproj.Proj(init='epsg:4326'), # source coordinate system
pyproj.Proj(init='epsg:32610')) # destination coordinate system
# g1 = gdf.loc[gdf['seg_index']==1000, 'geometry'].iloc[0]
# g2 = transform(project, g1) # apply projection
# create a dataframe with spectral values as column names , index value "seg_index" for merging geodataframe with Polygons
region_spectrals = []
for i in range(len(regions_red)):
seg_label = regions_red[i].label
g1 = gdf.loc[gdf['seg_index']==seg_label, 'geometry'].iloc[0]
g2 = transform(project, g1)
nir_zone_post = zonal_stats(g2, out_img_post_b08[0], affine=out_img_transform_post_b08, stats='mean', nodata=-999)
nir_zone_pre = zonal_stats(g2, out_img_pre_b08[0], affine=out_img_transform_pre_b08, stats='mean', nodata=-999)
dict_seg = {'seg_index': regions_red[i].label,
'red_value': regions_red[i].mean_intensity,
'blue_value': regions_blue[i].mean_intensity,
'green_value': regions_green[i].mean_intensity,
'nir_value': nir_zone_post[0]['mean'],
'red_value_pre': regions_red_pre[i].mean_intensity,
'blue_value_pre': regions_blue_pre[i].mean_intensity,
'green_value_pre': regions_green_pre[i].mean_intensity,
'nir_value_pre': nir_zone_pre[0]['mean'],
'area_m': regions_red[i].area * 0.31} #area in meters
region_spectrals.append(dict_seg)
df = pd.DataFrame(region_spectrals)
df.head()
# merge on seg_index
gdf2 = pd.merge(gdf, df, on='seg_index')
# #add in these two columns for labelling later
# gdf2['land_class'] = np.nan
# gdf2['burn_class'] = np.nan
gdf2.head()
def computeSI(b1, b2):
return (b1-b2)/(b1+b2)
def changedSI(SI_pre, SI_post):
return SI_pre - SI_post
# convert the keys to list
post_keys = ['blue_value', 'green_value', 'red_value', 'nir_value']
pre_keys = ['blue_value_pre', 'green_value_pre', 'red_value_pre', 'nir_value_pre']
perm = itertools.permutations(post_keys, 2)
perm_pre = itertools.permutations(pre_keys, 2)
post_SI_keys = []
for i, tup in enumerate(perm):
b1 = tup[0]
b2 = tup[1]
col_name = "SI_" + b1[0] + b2[0] + "_post"
post_SI_keys.append(col_name)
print(col_name)
gdf2[col_name] = computeSI(gdf2[b1], gdf2[b2])
gdf2.head()
pre_SI_keys = []
for i, tup in enumerate(perm_pre):
b1 = tup[0]
b2 = tup[1]
col_name = "SI_" + b1[0] + b2[0] + "_pre"
pre_SI_keys.append(col_name)
print(col_name)
gdf2[col_name] = computeSI(gdf2[b1], gdf2[b2])
gdf2.head()
compute dSIs
SI_keys = pre_SI_keys + post_SI_keys
SI_keys
band_combos = list(set([key.split('_')[1] for key in post_SI_keys]))
SI_combos = [ tuple([s for s in SI_keys if bcombostr in s]) for bcombostr in band_combos]
SI_combos
# now, add dSIs to dataframe
for i, tup in enumerate(SI_combos):
SI_post = tup[0]
SI_pre = tup[1]
col_name = "dSI_" + SI_post.split('_')[1]
gdf2[col_name] = changedSI(gdf2[SI_pre], gdf2[SI_post])
gdf2.head()
gdf2['land_class'] = np.nan
gdf2['burn_class'] = np.nan
gdf2.head()
gdf2.to_file("./data/segments.geojson", driver='GeoJSON')
fps_wv_post[5:6]
writeDatasets(fps_wv_post[5:6], fps_wv_post[5:6], fp_sent2_post, fp_sent2_pre)
#writeDatasets(fps_wv_post[8:], fps_wv_pre[8:], fp_sent2_post, fp_sent2_pre)